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ABSTRACT 


This  paper  presents  an  efficient  method  of  smoothing  steady-state, 
dissipative  hyperbolic  systems  with  one  spatial  dimension.  The 
observations  are  from  point  sensors  placed  on  the  system.  We 
show  that  under  realistic  stability  conditions  there  exists  a  family 
of  finite-dimensional  acausal  linear  systems  that  characterize  the 
frequency  domain  behavior  of  the  hyperbolic  system.  Using  this 
characterization,  we  develop  a  smoothing  algorithm  that  is  recur¬ 
sive  with  respect  to  the  sensors,  resulting  in  a  significant  decrease 
in  computational  complexity  relative  to  other  methods.  We  illus¬ 
trate  the  algorithm’s  performance  by  studying  the  smoothing 
problem  for  sound  waves  in  an  air-filled  pipe. 

1.  Introduction 

The  purpose  of  this  paper  is  to  derive  an  algorithm  for  the  linear  least- 
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squares  smoothed  estimate  of  inputs  or  state  variables  in  a  dissipative  hyper¬ 
bolic  system  described  by  a  vector  first-order  partial  differential  equation  with 
boundary  and  initial  conditions  [9]  .  Examples  of  such  systems  are  those 
described  by  the  telegrapher’s  equation,  the  vibrating  Timoshenko  beam  equa¬ 
tion,  and  the  wave  equation.  Our  smoothing  algorithm  can  be  used,  for  exam- 
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pie,  to  estimate  the  radiated  sound  from  a  vibrating  object  given  observations  at 
discrete  points  in  space  [3].  The  systems  we  consider  here  have  one  spatial 
dimension  and  are  operating  in  temporal  steady-state.  The  observations  are 
taken  by  Na  point  sensors  distributed  non-uniformly  across  the  system.  We 
assume  either  (a)that  the  observation  interval  is  long  enough  to  reliably  com¬ 
pute  the  Fourier  transform  with  respect  to  time,  and  to  cause  different  fre¬ 
quency  components  to  become  uncorrelated  with  each  other,  or  (b)  that  the 
relevant  random  processes  are  periodic  in  time  and  are  observed  over  an  inter¬ 
val  that  is  a  multiple  of  this  period. 

To  solve  this  smoothing  problem  we  first  Fourier  transform  the  observa¬ 
tions  with  respect  to  time.  We  then  have  a  set  of  uncoupled  spatial  smoothing 
problems  over  space,  indexed  by  the  frequency  variable  ,in  which  the  underly¬ 
ing  models  ore  finite-dimensional  well-posed  acausal  linear  systems  [4]-[5]. 
The  acausal  linear  system  smoothing  problems  are  solved  by  the  method  of 
complementary  models  [ l]  > [12]  after  which  one  may  use  an  inverse  Fourier 
transform  to  recover  the  estimates  as  functions  of  space  and  time.  The  result¬ 
ing  algorithm  is  recursive  with  respect  to  the  sensors,  and  thus  offers  a 
significant  decrease  in  complexity  relative  to  other  methods. 


2.  Dissipative  Hyperbolic  Systems 

In  many  signal  processing  problems,  one  has  measurements  of  the  output 
of  a  system  described  by  a  wave-like  (hyperbolic)  partial  differential  equation. 
Physically,  a  dissipative  hyperbolic  (DH)  system  [9]  is  a  model  for  a  wave  bear¬ 
ing  structure  that  has  internal  energy  loss  due  to  distributed  or  boundary 


damping.  Examples  of  such  systems  are  vibrating  strings,  beams,  transmission 
lines,  acoustical  and  electromagnetic  waveguides,  etc.  We  will  consider  DH 
systems  with  one  spatial  variable. 

A  DH  system  is  described  by  a  vector  first  order-system  of  partial 
differential  equations 

A  A 

-Jj m(x,t )  =  A{x)~m(x,t)+G{x)m(x,t)+e(x,t)  ,  x  €  [0,1]  ,  t>t0  (2.1) 
with  boundary  conditions 

Hom(0,t)  —  d^t)  ,  Hl  m(L  ,t)  =  d2(t)  (2.2) 

and  initial  condition 

m(x,t0 )  =  m0(x)  (2.3) 

where  m(x,t)  is  the  nXl  state  vector,  e(x,t)  is  the  input  field,  A(z)  is  asym¬ 
metric,  continuously  differentiable  matrix  with  constant  rank  r,  G(x)  is  a  con- 

0 

tinuous  matrix,  dL(t)  and  d2(t)  are  n/2Xl  boundary  inputs,  and  HQ  and  HL 
are  matrices  of  bounded,  linear,  causal,  shift-invariant  operators.  All  quantities 
in  Eqs.  (2.1)-(2.3)  are  real.  Note  that  according  to  the  Bochner- 
Chandrasekharan  theorem  [2],  the  boundary  operators  H0an.dHL  are  such 
that  their  operation  in  the  frequency  domain  is  multiplicative;  i.e.,  the  following 
transform  relations  hold: 

H0m(0,t)  HQ{ju)m(0,ju) 

Him(L,t)  Hi{ju)m(L  ,jw) 

where  and  Hi(jw)  are  complex  valued  n/2X  n  matrices. 


A  DH  system  will  satisfy  [9] 


m'(L,t)A(L)m(L,t)  -  m'(0,t)A(0)m(0,f)  <  0  ,  for  all*>f0  (2.5) 
when  dx  =  do  =  0.  These  conditions  ensure  that  when  e  =  dx  =  d2  =  0 

-|-l|m(x,t)|p  <  0  ,  for  allt>  t0  (2.6) 


where 


||m(a;,J)|p  =  /  m'(x,t)m(x ,t)dx 


To  see  this,  pre-multiply  Eq  (2.1)  (  with  e  =  0)  by  m'(x,t)  and  add  the  tran¬ 
spose  of  the  result,  to  obtain 

-£-m'(x,t)m(x,t)  =  ~-(m,(x,t)A(x)m(x,t))+m'(x,t)lG(x)+G,(x ) 


“  ^jA(i)]m(i,t) 


It  then  follows  that 


gjllm(x,t) |p  =  /  m'(x,t)(G(x)+G'(x)--£-A(x))m(x,t)dx 

+m'(L  ,t)A(L)m(L  ,t)~  m'(0,t)A(0)m(0,t)  <0 
An  example  of  aDH  system  is  the  damped  wave  equation 

vtt-  c2«„-H7u{  =  e 

with  boundary  conditions  (damped  supports) 

«<(0,0-*o«i(0,t)  =  dj(t) 

vt(L,t)+kLux(L,t)  =  d2(t),  k0,kL>0 
and  initial  conditions 

tt,(x,0)  =  gx{x)  ,ut(x,0)  =  g2(x ) 
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—  CU2  ( X ,  tj  j  7Uo(x,i)  — 


one  obtains 


dmx  0  c  q  mx 
~dt  m2  ~  cO  97  m, 


‘Iff3  lKL-pi 

2^0-7  m2  ^(lj 
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Note  that 


G+G'-  ~-A 
ox 
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and  when  dx  =  d2  = 


m\L  ,t)Am(L  ,t)~  m'(0,t)A.m(0,t)  —  -  2(k0mf  (0,t)+kL  mx  (L  ,t))  <  0. 

We  will  assume  that  the  DH  system  is  asymptotically  stable.  That  is,  if 
m(x,f)  is  the  solution  of  Eq  (2.1)  with  e  =  dx  =  d2  =  0  then  ||m(x,£)||  — *•  0  as 
t  — *  oo.  If  we  had  assumed  that  the  inequality  in  Eq  (2.4)  was  strict,  then  the 
inequality  in  Eq  (2.6)  would  be  strict  also,  giving  us  asymptotic  stability.  How¬ 
ever,  we  want  our  results  to  apply  to  systems  like  the  above  example  that  do 
not  have  a  strict  inequality  in  Eq  (2.4),  but  still  are  asymptotically  stable.  One 
should  note  that  a  system  that  has  normal  modes,  i.e.,  non-decaying  responses 
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to  initial  conditions,  is  not  asymptotically  stable.  In  practice,  however,  one 
always  has  dissipative  elements  in  the  system  and  these  elements  should  be 
retained  in  the  model  to  ensure  a  proper  formulation. 

Yet  another  stability  assumption  is  required  for  the  smoothing  problem 
when  the  A  matrix  is  singular.  In  this  case  it  is  shown  in  Appendix  A  that  each 
member  of  the  family  of  finite-dimensional  systems  comprised  of  those  state 
variables  associated  with  the  zero  eigenvalues  of  A  has  poles  with  only  nonposi¬ 
tive  real  parts.  In  order  for  the  results  of  this  paper  to  apply,  we  must  assume 
that  all  of  these  poles  in  fact  have  negative  real  parts.  Again,  if  the  inequality 
in  Eq  (2.4)  is  strict,  then  this  assumption  is  automatically  satisfied.  Physically, 
these  state  variables  correspond  to  damping  of  the  hyperbolic  system.  This  sta¬ 
bility  assumption  can  be  shown  to  ensure  that  an  acausal  linear  system 
representation  exists  and  is  well-posed.  The  well-posedness  of  the  acausal  linear 
system  then  implies  that  the  resulting  smoother  is  well-posed. 


3.  Problem  Statement  and  Construction  of  the  Acausal  System 

We  wish  to  determine  the  linear  least-squares  smoothed  estimates  of  the 
state  m(x,t)  and  inputs  e(x,t),dl(t),d2{t)  of  the  DH  system  (2.1)-(2.3)  given 


observations 


?*(*)  =  Cm(xklt)+wk(t) 


at  specific  points  xk  along  the  system,  where  C  is  p  X  n  and 


t  e  [-r/ 2  ,  772],  k  =  1,2, 


0<x1<x2<  '  ’  ’  <xNt<L 
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We  assume  that  e(x,t)ltf1(t),J2(0  and  observation  noise  «;*(£)  are  zero  mean 
and  wide-sense  stationary  in  time,  that  t0  =  -oo  and  m0(x)  =  0.  The  signals 
m(x,t)  and  yk ( t)  are  then  also  wide-sense  stationary  .  We  also  assume  that 

Ee(x ,t)e'(z,s)  =  Q(x,t-  s)S(x-  z) 

Ewi(t)w/(s)  =  Rw(t-s)S,y 
Ev(t)v'(s )  —  n„(f-s)  where  v(t)  =  [<//(<)  d^t)]’ 


Ee(x  ,t)wi'(s)  —  Ee(x,t)v,(s)  =  Ew^fyv^s)  =  0 


We  will  use  a  Fourier  series  expansion  (in  time)  of  the  signals  yk{t), 

m(x,t),e tc.,  over  the  interval  [- T / 2  ,  r/2],  and  denote  the  Fourier  coefficients 

by  yk(ju}),m(x ,jw),  etc.,  where  u  —  2vl/T  ;  l  —  0,±1,±2,  •  •  •  .  These 

# 

coefficients  can  be  computed  using: 

r/2 


yk{j(j)  =  -=  J  yk(t)e->“‘dt 
1  -  r/2 


Note  that  this  integral  can  be  evaluated  at  u  —  2 irl/T  using  discrete-time  data 
and  an  FFT  if  the  signal  is  band-limited.  The  stationarity  assumption  implies 
that  the  Fourier  coefficients  are  uncorrelated  at  different  frequencies  provided 
one  of  the  following  conditions  holds  [8]  :  (1)  the  covariances  Q(x,r),  Rw(r), 
n„(r )  are  periodic  in  r  with  period  T  ;  or  (2)  the  observation  interval  T  is  long 
and  the  covariances  go  to  zero  as  r  — ►  oo.  We  note  that  in  many  vibration  and 
acoustical  problems,  the  signals  (and  hence  their  covariances)  are  periodic.  In 
the  following  we  will  assume  that  one  of  these  conditions  holds,  in  which  case 
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coefficients  of  Q(i,r),nt)(r),  Rw{r),  which  when  T  is  large  can  be  approxi¬ 
mated  by  dividing  the  spectral  densities  of  e  ,v,  and  w *  by  T.  We  will  assume 
that  Rw(ju>)  is  invertible. 


If  A  (a:)  is  invertible  for  all  x,  we  can  write  Eq  (3.1)  as 

Q 

——m(x,ju)  =  A(x,ju>)m(x,ju)+B(x)e(x,ju>)  (3.4) 

ox 

where 

A(x,juj)  =  A~l{x){ju}I-  G(x)} 

B(x)  —  -  A_1(x) 

Eqs.  (3.4), (3.2), (3.3)  are  a  family  of  acausal  linear  systems  indexed  by  u>.  It  is 
shown  in  Appendix  A  that  a  similar  acausal  linear  system  representation  can  be 
obtained  even  when  A(x)  is  not  invertible,  provided  the  stability  assumptions 
discussed  earlier  hold. 

We  must  consider  the  well-posedness  of  the  acausal  linear  system  descrip¬ 
tion  of  the  DH  system.  An  acausal  linear  system  is  well-posed  if  there  exist  no 
nonzero  solutions  to  an  undriven  system. 

Theorem :  The  acausal  linear  system  of  Eqs  (3.2)-(3.4)  is  well-posed  for  all  cj  . 

Proof :  Suppose  there  exists  an  w0  such  that  Eqs  (3.2)-(3.4)  are  not  well- 
posed.  There  then  exists  m(x satisfying 

Q 

— -m{x,jw0)  =  A{x,ju0)m(x,ju0) 

OX 

VQm(0,j<JJ0)  +  VLm(L,ju0)  =0 

or,  equivalently, 

Q 

juj0m(x  ,ju>Q)  =  A(x)~—m(x,ju>0)+G(x)m(x,juj0) 


#o(M))m(0>iw0)  *=HL(ju0)m(L,j<jj0)  =  0 


It  is  easily  checked  that  if  tf(x,t)  =  eJUotm(x, ju>0),  then  *(x,t)  satisfies  Eqs 


(2.1)  -  (2.2)  with  zero  inputs: 


■j^<i/(x,t)  =  A(x)j-V(x,t)  +  G(x)V(x,t) 


H0*(0,t)  =HL*{L,t)  =0 


Since  ||'I'(x,i)||  does  not  go  to  zero  as  t—>o o,  we  have  a  contradiction  of  the 


asymptotic  stability  assumption. 


4.  Smoothing  the  Acausal  Linear  System 


In  this  section  we  show  how  to  solve  the  smoothing  problem  for  our 


acausal  linear  system.  Although  this  paper  concentrates  on  DH  systems,  many 


parabolic  type  equations  can  also  be  written  in  an  acausal  linear  system  form. 


The  smoother  is  derived  by  means  of  complementary  models,  introduced  by 


Weinert  and  Desai  [12],  and  extended  by  Adams  ,  Willsky  and  Levy  [ l] .  The 


derivation  differs  significantly  from  that  in  [l]  due  to  the  possible  singularity  of 


Q  and  IIU  ,  and  the  fact  that  the  observations  are  discrete.  In  what  follows,  the 


w  dependence  in  Eqs.  (3.2)-(3.4)  will  be  suppressed. 


A  solution  to  Eqs.  (3.2), (3.4)  is 


m(x)  =  $(x,0)F~  lv+f  G(x  ,z)B(z)e(z)dz 


(4.1a) 


where  the  state  transition  matrix  $  satisfies 


— $(*,*)  =  A(x)$(x,z),  $(z,z)  =  I 


(4.1b) 


and  the  Green  function  is  given  by 


G(x,z ) 


#(x,0)F-17o$(0,a)  if  2<x 
-*(x,0 )F~1Vl*(L,z)  \tz>x 


(4.1c) 


and  the  matrix  F  satisfies 


J®*-  V0+Vl*(L, 0)  (4. Id) 

The  well-posedness  of  the  foregoing  acausal  linear  system  guarantees  the  inver- 
tibility  of  F  [4]-[5] . 

If  rank  IIV  =  g x  we  can  write  a  full  rank  factorization  of  nv  as 

Uv  =  MM * 

where  M  is  »X  g x;  in  other  words, 

v  —  Mfi  ,  Ejxm  *  =  Ju 
Similarly,  if  rank  Q{x)  =  g2,  we  can  write 

Q(x)  =  S(x)S*(x) 

where  S(x)  is  nXg2;  thus 

e(x)  =  S{x)p(x)  ,  Ep(x)p*(z)  =  IJaS{x-z) 


Using  Eq  (4.1a),  we  can  write  Eq  (3.3)  as 

L 

yj  =  C$(xj,0)F~lM M+f  CG(.  z)B(z)S(z)p(z)dz+Wj  (4.2) 

o 

Eq  (4.2)  relates  the  observations  to  the  underlying  variables  m, p(  ), {toy},  which 
together  span  a  Hilbert  space  H.  If  Y  is  the  Hilbert  space  spanned  by  the 
observations  {y;},  then  it  is  shown  in  Appendix  B  that  yc(  •)  and  6  defined 


below  span  the  orthogonal  complement  Yk 


Vc(x)  ~  p(x)-S*(x)B*(x)\(x)  (4.3a) 

0  =  /i-M*{F*}-1(X(0)-$*(L,0)X(L))  (4.3b) 

where 

K 

X(z)  =  2  G*(xk,x)C*Rw~1wk  (4.4) 

£  =1 

Now  Eq.  (4.4)  implies 

-^-X(ar)  =  -A*(x)X(x)  ;  xj&Xj  (4.5a) 

\{xr)  =  \(xj+)+C*Rw~1Wj  ,  j  —  1,2,  •  •  •  ,Na  (4.5b) 


In  order  to  make  Eqs.  (4.5)  equivalent  to  Eq.  (4.4)  a  boundary  condition  of  the 
following  form  is  needed: 

» 

Kq\{0)+Kl\{L)  =  0  (4.5c) 

where  K0  and  KL  are  nXn  matrices  .  Eqs.  (4.1c)  and  (4.4)  imply  that  Eq. 

(4.5c)  holds  if 


I(0V*  -  ICL  Vl-  0  .  (4.6) 

Furthermore,  if  K0$  *(L  ,Q)+KL  is  invertible,  then  Eqs.  (4.5)  will  be  well- 


posed.  If  we  take 


*o-  wr* 


Kl  =I-VZ{F*}-'**(L,0) 

then  both  the  invertibility  condition  and  Eq  (4.6)  will  be  satisfied.  With  this 
choice  of  K0  and  K i,  the  acausal  linear  system  (4.3),  (4.5)  is  complementary 
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to  (3.2)-(3.4). 

Solving  Eq.  (4.3a)  for  p(x)  and  substituting  into  Eq.  (3.4)  gives 
d 


(4.8) 


— m(x)  =  A(x)m(x)+B(x)S(x)(S  V)£'(*)M*)+ye(*))  (4.7) 

Likewise,  solving  Eq.  (3.3)  for  wj  and  substituting  into  Eq.  (4.5b)  gives 
\(xr)  =  \(xJ+)+C*Rw~1(yr  CmiXj)) 

Multiplying  Eq.  (4.3b)  by  M  and  re-arranging,  we  get 

v  =  M0+nv{F*}-  1(X(0)-  $  *{L  ,0 )X(L )) 

Combining  Eqs.  (3.2),(4.5),(4.7)-(4.9)  gives  the  Hamiltonian  system 


(4.9) 


d  m(x, 
dx  X( x  , 


\A  BQB' 
0  -A* 


[>(  *?  ’  1  ^Xj  ^ 4-10a^ 


X(xr)  _  X(x;+)-  C*Rw~lCm(xj)+C*Rw-1yj 


-i. 


(4.10b) 


K0  -n v{F*rl 
0  K0 


MOR 

U(0)  f 


VL.nv{F'}-l*\L,0) 

0  kl 


\m 


(4.10c) 


By  projecting  the  random  quantities  in  the  Hamiltonian  onto  the  subspace  Y 
spanned  by  the  observations,  we  obtain  the  estimate  Hamiltonian: 


d 

m(x) 

A  BQB*' 

m(x) 

dx 

X(x) 

0  -A* 

X(x) 

,  Xf^Xj 

\(xr)  =  X(xy+)~  C*Rw~1Crh{Xj)+C*Rw~  lyj 


(4.11a) 

(4.11b) 


0  = 


v0  -nv{F*yl 
o  icQ 


o 

■  Jj 

1 

4- 

jlX(O)  j 

I 

vL  nP{F*r^*(L,o) 

o  kl< 


m(L ) 

\{L) 


(4.11c) 
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where  ‘denotes  projection  onto  Y.  To  obtain  e(x)  and  v,  project  onto  Y  in  Eq. 
(4.3a)  (after  multiplying  by  S),  and  in  Eq.  (4.9): 

e(x )  —  Q(x)B*(x)\(x)  (4.12a) 

v  =  nv{F*Yl(\(0)-$*(L  ,0)\(L))  (4.12b) 

To  prove  that  Eqs  (4.11)  are  well-posed,  assume  that  y}-  —  0  for  all  j.  Now 
fh(x)  is  the  linear-least  squares  estimate  of  m(x)  based  on  observations  that 
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are  all  zero,  so  m(x )  —  Em(x)  =  0,  the  last  equality  following  from  the  well- 
posedness  of  Eqs.  (3.2)-(3.4).  Moreover,  since  X(x)  satisfies  the  unforced  ver¬ 
sion  of  Eqs.  (4.5),  which  are  well-posed,  X(x)  =0  for  all  x  ,  and  hence  Eqs. 
(4.11)  are  well-posed. 

5.  Recursive  Solution  to  the  Estimate  Hamiltonian 

To  solve  the  estimate  Hamiltonian  (4.11),  we  will  first  diagonalize  the 
dynamics  with  the  following  change  of  variables  : 
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=  T(x) 
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P(x)  -(«,(*)+«, (*))-• 
and  9 f{x)  and  9b(x)  satisfy 


_  ~  ^ fA+A*0 f+8 jBQB*0 j 


— —  8bA+A  6b-0bBQB*8b  ,x^x/ 
with  any  positive-definite  boundary  values  6/(0)  and  6b(L).  The  new  variables 
^  /  and  ^  t  satisfy 

_d  f*/(x)l  f*/  01|J'/(I)1  , 

&[*((*)]  [o  Atht  tlx)\<  X^XS  (5-la) 

\v?  n°i  +!v/  vf]  i[\LL\  =0  (s.ib) 


where 


A/(x)=-(A\x)+9/(x)B(x)Q(x)B\x)) 
Ab(x)  =-(A\x)-6b(x)B(x)Q(x)B*(x)) 


f  o  o  _  f*ro+n,{/,V'*(o)  K0-nv{n-^/(o)]rF(o)  o 

\vf  v»]  -  -Ko0b(O)  Ko0/( 0)  0  P(O) 


if]  = 


VL-nv{Fy'*\L,0)6b(L)  Vl+Uv{F*}-^\L,0)6/(L)\P{L)  0 


-KJ»(l) 


KJf(L) 


0  P(L) 


We  must  now  specify  the  evolution  of  'i ’ /,<!/ b,9/,8b  at 

*  *  X/  ,j  =  1,  •  •  •  ,Nt.  If  we  choose 

6f(xi+)^9f(xr)+C*RV)-lC 


Ob(xr)=*0b(xj+)+C'Rw-'C 


then  one  can  show  that 


*/(*;+)  -‘*fi*r)+C*Rv~l9j 


(5.2a) 


=  'if  i(x]+)+C*Rw~1yj  (5.2b) 

Eqs  (5.l)-(5.2)  are  in  a  form  that  can  be  solved  recursively.  In  terms  of 
y (0)  and  ^  t,(L)  ,  a  solution  to  Eqs.  (5.1)-(5.2)  is 

♦/(*)  -*4f(*.0)*/(0)+*7(*)  (5.3a) 

*»(*)  *  ^(x,I)^6(I)+^60(®)  (5.3b) 

where 

-^r*Af(x’ o)  =  ;  ^(M)  =  / 

-&*At(x’L)  =  aU*)*a>(x>l)  ;  *aAl>l)  —  i 


-£-*/(*)  =  A7(*)«  J{x)  ,  x^Xj  (5.4a) 

-£-*»°(*)  =  A6(x)¥?(x)  ,  x^x;-  (5.4b) 

*7(0)  =  *?(£)  =0  (5.4c) 

*?(*/+)  =  #?(*,- HCX'1*  »  >  -  1.  *  •  •  .*•  (5.4d) 

#?(*/-)  =  ♦ftxy+J+CX"1^-  ,  y  =  1,  ■  •  •  ,iV,  (5.4e) 


Setting  x  —  L  in  Eq.  (5.3a)  and  x  =  0  in  Eq.  (5.3b),  and  using  Eq.  (5.1b)  we 
obtain 


*/(0) 

*»(£) 


-VOW^I+Wt0)) 


t'*r'  ^vp*  — 


Fn  =  [vZ+Kf-^cz , o) :  vf+v?*A(o,i)] 


A  recursive  solution  to  Eqs  (5.1)-(5.2)  is  therefore  given  by: 


0  $Ai(x,l) 


'*/“(*) 


*/(*) 

*»(*) 


^|*4°(x)J 

in  conjunction  with  Eqs.  (5.4). 

We  now  show  that  F ^  is  invertible.  When  yj  —  0  ,  j  —  1,  •  •  •  ,Nt,  anj 
solution  to  Eqs  (5.1)-(5.2)  will  satisfy 

*/(*)  =  *A,(x>  0)01 


where 


*»(*)  =  *aSx’L)P  2 


F  —  0 

Ffl  02  ~  ° 


If  Fjt  were  not  invertible,  we  could  generate  non-zero  solutions  to  Eqs.  (5.1)- 
(5.2)  and  hence  to  Eqs  (4.11)  (via  T^x))  with  yj  —  0,  contradicting  the 
well-posedness  of  Eqs  (4.11). 


6.  Smoothing  Error  Covariance 

We  wish  to  calculate  the  time  domain  error  covariance 

Pm{xA  —  E[m{x,t)m\x,t)\ 

where  m  =  m-  m.  Using  Fourier  transforms,  one  can  show  that 


where 


Pm(x,t)=Pm(x, 0)  =  £  Pm(x,j2nl/T ) 

l  =  _  00 


PmiX’ju)  —  E[rh{x,juj)m\x,ju)\ 


A  dynamical  equation  governing  m(x,ju))  is  obtained  by  projecting  the  Hamil¬ 
tonian  system  (4.10)  onto  Y^,  and  using  the  definitions  of  yc  and  6: 

d  ’m(x)  ]_  P1  lfm(ar)  ],  pel  , 

~dx  J- X(x)  [o  -A*  J|-X(4f"Lo  J’  ***/  (6>la) 

-Z{xr)  =  -1^+)-  C*Rw-lCm{Xj)-  C'Rv-'wj  (6.1b) 


\v]_  V°  Uv{F  }~  fm(0)  l 

jo  J  0  IC0  |_X(0)J4' 


_  ^.-n,n,|L(0)l  \Vl  nv{F*yl**(L>0)  Ir^j 


-X(L) 


(6.1c) 


Here  again  the  w  dependence  has  been  suppressed.  Note  that  these  equations 
are  similar  to  Eqs.  (4.11)  and  can  therefore  be  solved  using  the  same  diagonal¬ 


izing  change  of  variables.  Thus,  if 
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eb{x) 


-  Ax) 


\m(x)  ] 


where 


e/(x) 


eP(x) 


C(X)  =  */»(*)  Vi*-  vh/°(L)-  nvwi*  eyix) 


-j-e?{x)  =  Af(x)e?{x)+9f{x)B(x)e{x)  ,  x^Xj- 


-4-e^x)  =  Ai(x)ei°(x)+^(X)B(2 :)€(x)  *  x7^xj 


«/(0)  =  e,°(£)  =  0 


«/(^y+)  -  )-  C*Rw~1vij  ,  y  _  i, . . .  jA, 

)  -  «.“(*/+)-  eV,j ,  j  =  i,  .  .  .  jJV> 


*/?(*)  = 


Letting 


*^(*,0J  0 

0  **»(*,£) 


e(l)-^w)<;(l)e*'(l)1“p 


|®u(*)  eX2(z) 
e2i(*)  e22(x) 


pm{x,ju>)  «  ^(^)(e11(x)+e22(I)+©12(a.)+©2i(a.)jP(i) 
With  the  following  definitions: 

n7(*)  »£[€/(*)«/•(*)} 

n*(*)  ~  E[t?{x)e?\x)\ 
RMx,y)=E{ef°(z)e<>\y)l 

0(x)  can  be  written  as 

e(l)  “  */»(*) V fi.+»fn,(£ ) vf+vfRtfi  ,0) k,“-+ v.-'n.fo) n°  + 

n7(x)  o  * 

o  n*(») 
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J(x)  =  F/^/(L,x)n/(x);i2/6(L,i)]fn°[^(a:,0);n6(x)<J>;t(0,x)] 

By  direct  evaluation,  one  can  obtain  the  following  formulas  for  R ^  ,  TIf  ,11$: 

X 

R/t(x,v)  —  -/ $A,(x><r)ff/(a)B(‘7)Q(a)B*(ar)ffi(a)*At(y ^)da,  x>y 
y 


-~Uf(x)  =  AjTlf(x)  +  Tlf(x)Af  +  OfBQB*df ,  sj^afy 

ax 

4-^b{x)  =  AbUb(x)  +Ub(x)Ab*~  6bBQB*9b  ,  x^Xj 
ax 


n,(o)-nb(L{)-o 

nf{Xj+)  =n;(ir)+cVlc  ,  j  -  i,  ••  •  ,jv, 

II4  (ly-  )  =  II 4  (ly+J  +  C1  i?*,-  ,  j  —  1,  >-Nj 
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7.  Separable  Boundary  Conditions 

In  this  section,  we  show  that  if  the  boundary  inputs  db(t)  and  d2{t)  of  the 
DH  system  are  uncorrelated,  and  if  the  resulting  block  diagonal  II  „  is  inverti¬ 
ble,  the  smoother  and  error  covariance  formulas  simplify  considerably.  Under 
these  conditions,  V^II"1!^  =  0  ,  in  which  case  the  acausal  linear  system  has 
separable  boundary  conditions.  To  examine  the  filter  and  smoothing  error 
covariance  for  separable  boundary  conditions,  premultiply  Eq.  (4.11c)  by 

VoX'1  -*'(L,0)' 

VlU;1  i 

giving 


r0Xln  -r 

m(0) 

4- 

0  o' 

’m(L) 

0  0 
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yi n.-'n  / 
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If  we  choose 

W°)  =  v*n;'v0,eb(L)  =  v[n-v'vL 

then  the  change  of  variables  using  T(x)  produces  the  boundary  conditions 


so  that  ¥  /  — 


*,(0)  =  *&(L)=0 

tyf  and  'I'  b  —  ^b.  Therefore 

m(x)  =  P(x)'i'J(x)+P(x)'i!  ?(ar) 


\(x)  =  0f(x)P(x)V¥(z)-$i(x)P(x)<!rf(x) 

Furthermore, by  proceeding  in  the  same  fashion  with  Eq.  (6.1c),  we  get 


so  that 


It  follows  that 


e,( 0)  -  V*QU~vlv  ,eb(L)  =  itfl;1* 

ef(x)  mm  $Al(x,0)VZn;lV  +  C?{x) 
eb  (*)  =  *Ah(xiL)VtTl;xv+c?{x) 

012(i)  =  ©21(2:)  =  0 


©11(1)  =  0)tf/(0)^^(x,0)+n/(x)  =  6f{x) 

en(x)=*Ai(*,L)6i(L)<l>l(x,L)+  Ylb(x)  =  6b(x) 

Therefore,  Pm(x,jw)  =  P(x). 

8.  Algorithm  Complexity 

We  will  examine  the  computational  complexity  of  the  algorithm  presented 
in  Section  5.  For  computational  purposes,  we  assume  that  the  interval  [0,L]  is 
partitioned  into  X  regions  ,  and  that  an  FFT  with  Cl  frequencies  is  used.  Table 
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1  1  gives  the  complexities  of  the  major  steps  necessary  to  compute  the  estimate 

* 

!  »»(  v)  . 


TABLE  1  Complexity  of  the  Smoothing  Algorithm 


Step  Complexity 

Fourier  transform  the 
observations 

0(Q  NaP logQ  ) 

Compute  * 

0(  Q  [ Xn3+N„  np  +N,  n2] ) 

Com  pute  m  ( x ,  ju> ) 

0(QXn2) 

Inverse  transform 

ih(x,joj) 

0(f2  nXlogfl  ) 

As  can  be  seen  from  Table  1,  the  overall  complexity  of  the  algorithm  is 


s.' 

V 

a 

W 

w" 


v; 


0(11  [X(n3+nlogfl  )4-X4(np+7J2+plogfJ  )]) 

For  comparison,  if  a  Wiener  smoother  is  used 

m(x,t)  =  ^  1  {SmY(x,ju)Sp^(joj)  Y(ju)} 

where 

Y{ju)  =  [vf(ju),  ■  •  •  ,yfc(ju))T 

SmY(x,ju)  is  the  cross-spectral  density  between  m(x,ju>)  and  Y(ju),  and 
SyyUu)  "1S  t^ie  spectral  density  of  Y(ju );  then  the  complexity  of  the  recon¬ 
struction  is 

0( n  [Ntp(nX+N-p2+ logfl  )+«Xlogf)  ]) 

An  example  using  this  type  of  approach  can  be  found  in  [11].  We  see  that  the 

algorithm  presented  in  this  paper  is  most  advantageous  when  the  number  of 

sensors  are  large  (linear  in  Ns  versus  cubic  in  Na). 
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9.  Example:  Sound  Waves  in  an  Air  Filled  Pipe 

As  an  example  of  the  use  of  our  algorithm,  we  consider  the  problem  of 
estimating  the  particle  velocity  and  sound  pressure  levels  inside  an  air  filled 
pipe  with  rigid  walls,  given  observations  of  sound  pressure  levels  at  discrete 
points  along  the  pipe.  We  are  interested  in  estimating  the  velocities  and  pres¬ 
sures  in  the  section  of  the  pipe  from  x  —  0  to  x  =  3  meters,  and  we  assume 
that  forward  traveling  waves  enter  from  the  x  =  0  boundary,  and  backward 
traveling  waves  enter  from  the  x  =  3  boundary.  The  particle  displacement  is 
assumed  to  satisfy 

■4>tt{x,t)  =  cfy„(x,t)+e(a:,0  ,  *  6  (0,3) 


ipt(0,t)-  cipx(0,t )  =  <^(0 


^t(3,0+«^*(3,0  “  rfa(0 

where  c  is  332  meters/second  ,  e  is  a  noise  term  accounting  for  yielding  of  the 
pipe  walls,  and  dx  and  do  are  waves  entering  the  pipe  section.  The  observations 
are 


Vj(  0  =  -  PoC’tlfx(xjtt)+Wj(t)  ,  j  =  1,  •  •  •  ,  Ns 
where  p0  is  the  density  of  air  (1.29  kilogram /meter  3),  and  Wj(t)  is  observation 

noise.  In  DH  form  this  system  becomes 
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n»i(i,i) 

o  c' 

d 

m  x(x,t) 

'  0  1 

dt 

m2(x,t) 

. 

c  0 

. 

dx 

m2(x,t) 

+ 

e(i,£)  J 

Mil! 


m2(0,t) 


-  dx{t)  ,  [1:1] 


mi(3,£) 

m2(3,£) 


d2(t) 
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Vj{t)  -  [-p0c  O}  „,2(xj|f)  jWO  J  =  1.  •  •  •  >Ns 

where 

m\(xyt)  =  cijjx{x,t) 
m2(x,t)  =  i/>t(x,t) 

We  make  the  following  statistical  assumptions 

R~1(ju)  —  10“ 2  pascals2 

|7  °1 

n,(jw)  =  10“ 6  q  ^  meters2 j seconds2 


[0  °1 

Q{jw )  =  10“ 2  q  j  meters2 /seconds* 

for  lo  =  2it  l  radians  /second,-,  l  =  ±1,±2,  •  •  •  ,±500.  Figures  1  and  2  show 
the  smoothing  error  covariance  for  the  sound  pressure  and  particle  velocity 
respectively,  as  a  function  of  the  number  of  sensors  {Ns)  uniformly  distributed 
along  the  pipe.  Figures  3  and  4  show  the  smoothing  error  covariances  as  a 
function  of  frequency  for  a  pipe  with  5  sensors.  In  these  figures,  the  x-axis  is 
frequency  (from  0  to  1000  it  radians  /second)  and  the  y-axis  is  length  along  the 
pipe  from  0  to  3  meters.  One  can  see  the  effects  of  spatially  sampling  the  sound 
field.  In  Figure  4,  the  error  covariance  maximum  occurs  at  w  =  664rr  radians 
/second.  Since  the  spatial  sampling  frequency  for  5  sensors  is  4a-  radians  / 

meter  the  error  covariance  maximum  occurs  when  the  wavenumber  k  =  —  of 

c 

the  sound  waves  matches  the  Nyquist  sampling  frequency.  Figure  5  shows  the 
actual  and  reconstructed  time  waveforms  for  the  sound  pressure  level,  using  8 


sensors.  The  actual  and  reconstructed  pressure  field  as  a  function  of  space  at 
the  frequency  u;  =  500jt  radians  /second  appears  in  Figure  6.  Figure  7  displays 
the  actual  and  reconstructed  particle  velocity  as  a  function  of  space  at  the  same 
frequency. 

10.  Concluding  Remarks 

The  input  estimates  can  be  interpreted  as  the  result  of  a  generalized  Bom 
inversion  procedure.  For  instance,  in  the  case  of  the  ID  wave  equation, 
e(x,ju ;)  will  be  the  Born  approximation  to  the  wave  speed  variations  in  an 
inverse  scattering  experiment  [7].  In  such  problems,  one  may  update  the  wave 
speed  function  in  an  iterative  fashion.  The  approach  used  in  this  paper  to 
derive  the  smoothing  algorithms  is  based  on  using  a  frequency  domain  two 
point  boundary  value  problem  to  describe  the  system’s  dynamics.  A  related 
approach  to  characterizing  a  vibrating  system’s  dynamics  is  given  in  [6],  where 
variations  in  the  system’s  parameters  are  assumed  to  occur  at  discrete  points 
along  it’s  length,  giving  rise  to  a  constant  diagonal  A  matrix  .  This  type  of 
model  can  be  handled  using  the  algorithms  developed  in  this  paper.  In  both  of 
these  approaches,  one  can  interpret  the  boundary  conditions  at  the  endpoints  of 
the  system  as  describing  the  reflection  and  transmission  coefficients  of  the 
hyperbolic  system.  In  the  DH  case,  the  reflection  and  transmission  coefficients 
arise  in  a  natural  way  when  the  A  matrix  is  diagonal.  Note  that  for  many  dissi¬ 
pative  systems,  there  does  not  exist  a  discrete  set  of  spatial  eigenfunctions,  so 
that  a  modal  expansion  of  the  dynamics  and  observations,  a  technique  used 
quite  often  in  distributed  parameter  filtering  and  control,  is  not  in  general  appli- 
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cable  to  the  systems  discussed  in  this  paper.  The  frequency  domain  description 
of  hyperbolic  systems  with  2  and  3  spatial  dimensions  involves  distributed 
parameter  acausal  linear  systems.  Efficient  smoothing  algorithms  for  the  2-D 
wave  equation,  for  example,  can  be  developed  in  this  way  [10]. 
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Appendix  A 

To  study  the  case  when  A(ar)  is  not  invertible,  we  put  Eqs  (2.l)-(2.3)  into 
a  canonical  form  which  separates  the  distributed  parameter  states  from  the  so- 
called  local  states.  The  local  state  variables  typically  correspond  to  damping 
forces  acting  on  the  distributed  parameter  system.  These'  damping  forces  may 
be  due  to  external  inputs,  corresponding  to  an  active  control  system,  or  the 
forces  may  be  passive,  such  as  structural  damping  in  a  beam.  Phillips  [9] 
proves  the  existence  of  a  family  of  orthogonal  matrices  { U(x )  ;  0<x<L}  with 
absolutely  continuous  elements  having  square  integrable  derivatives  such  that 


U\x)\(x)U(x)  =  [° 

•  < 


where  positive  definite  and  r  by  r.  Phillips  gives  an  explicit  algorithm 
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for  calculating  U(x).  With  the  change  of  variables 


y{x,t)  =  U'(x)m(x,t ) 


Eqs.  (2.l)-(2.3)  become 


-  di(‘)  •  UiL)„(L,t)  =  d2{t) 


y(x,t0)  =  U\x)mQ[x) 

The  stability  assumptions  are  unchanged,  because 

U'\V,+U',AV+U'G'V+V'GU-£°0  Am°(i)  -  t/'(G+G'--^A)V  <  0 

0  A22°(L)  =  m'(L,t)A(I)m(L,0 

y'(0,i)  0  Ao!(0)  »(°»0  =  m,(0,OA(0)m(0,f) 

- 

We  therefore  assume  that  Eqs  (2.1)-(2.3)  have  the  following  canonical  form 


d  _  0  0  <9  <"'12 

dt  m<z(x,t)  0  A22(i)  nio(x,t)  G  21  ^*22  |Jm2(i,0 


ej(a;,0 

+  e2(a;,0 


(A.l) 


[^01  :  tfoa]  micolo  =  d ’  K1  :  m2(l!o  ~  ^(0  (A-2) 


^  ;• . 

'  S*  ->* 
0  *  ■  • 


/..•  V 

w.v 

aW 


V  V% 

<v  \ 
Jrf- A- 


••  v  v 

rj 
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m0{x)  =  0 


Fourier  expanding  Eqs  (A.1),(A.2)  at  w  =  2nl/T,  l  =  0,±  1,  •  •  •  ,  we  get 

0  0  5  m|(i,jw)  _  juI-Gn  -G 12  mx(x,ju}) 

0  A22(x)  m2(x,joj)  —  G2X  jwl-  G22  m2(x,j<j) 


e^xjw) 

e2(x,j'w) 


(A. 3  a) 


j/f01(ia;)  ;  if02(jw)  m^o,jw)  —  (A.3b) 

-<><*>)  (a.3c) 


The  stability  assumption  (2.4)  implies  that 


<?n+<?u'  <  0 


(A.4) 


The  first  row  of  Eq  (A.3a)  gives 


(juiI-Gn)mx{xtj<jj)  —  G12m2(x,ju)+ex(x,jw) 


(A.5) 


From  Eq  (A.4)  we  see  that  the  eigenvalues  of  Gn  have  negative  or  zero  real 
parts.  The  stability  assumptions  that  were  discussed  in  Section  2  eliminate  the 
possibility  that  these  eigenvalues  are  on  the  imaginary  axis.  Solving  Eq.  (A.5) 
for  mx(x,ju> )  and  substituting  into  Eqs.  (A.3)  gives 


—m2{x,jw)  =  A 22  [jwJ-  G22-  G2iUuI~  Gli)~lGlJt\m2{x,ju) 
-  A221(Gf21(yw/- Gn)-le1+e2) 


(A.6) 


[Hol{jwI-Gn{0))-lGx2{Q)+HQ2)m2{0,jw)  =  dx(ju)  (A.7a) 
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[HqlUuI-  g\\{L  ))~lG\i{L  )+HL2]mi(L  ,ju)  =  d2(joj)  (A.7b) 

In  deriving  Eqs.  (A.7)  we  have  assumed  that  e(0,/o;)  =  e(L  ,ju>)  —  0.  These 
equations  are  now  in  the  form  of  Eqs  (3.2), (3.4). 


Appendix  B 

In  this  appendix  we  will  show  that  ye(-)  and  0  defined  in  Eqs  (4.3)  span  Y J- 
,  the  subspace  of  random  variables  orthogonal  to  Y,  so  that  Y<9Y^-  is  a  direct 
sum  decomposition  of  H,  the  underlying  Hilbert  space  generated  by 
{ p  ,  p(x )  ,  0 <x<L  ,  Wj  ,  j  —  1,  •  •  •  Na}.  We  introduce  rjj  and  z j  defined  by 

Vj  —  R»x,2u>j  ,  */  —  Rw  1/2Vj 

so  Eq.  (4.2)  can  be  rewritten  in  an  obvious  operator  notation  as 


z  —  Ffi+Gp+rj 


(B.1) 


where 


z  =  •  •  •  zNy  i\  =  •  •  •  nN.r 

If  a  €  H ,  its  projection  onto  Y  is  denoted  d  and  its  projection  onto  Y^-  is  a. 
Decomposing  (p  ,  p  ,  tj)  gives 


where 


Therefore, 


If  we  define 


P  =  p+p  ,  p  *=  p+p  ,  rj  =  ff+rj 

r\  —  R~lz  ,  p  ~  G*fj  ,  p  =  F*rj 
R ,  =  [FF'+GG'+/j 


P  S  p-G*fi+G*rj  ,  p  =  h-F*v+F*V 


ye  —  P-  G *n  ,  6  —  p-  F*i j 


(B.2) 


R*  .*•  <•  *■? ..  •  '  y  . ■ «  w.  m  v. * 

MHJGtfCiOtXT 


t*  i 


yc  —  p—  G  rj  =  (I+G*G)p+G*Ffi 
6  =  fi-F'rj  =  F*Gp+{I+F*F)p 
where  we  have  used  the  fact  (  see  Eq.  (B.l)  ) 


(B.3a) 

(B.3b) 


0  =  Ffi+Gp+ij  (B.4) 

It  is  clear  from  Eqs.  (B.3)  -  (B.4)  that  (ye  ,  6)  uniquely  determine  {p,  ,  p  ,  ij) 

thus  (y  ,  yc  ,  0)  uniquely  determine  (p  ,  p  ,  ij).  Note  also  that  Eq.  (B.3) 

implies  ye  €  and  6  €  YK  To  verify  Eqs.  (4.3)-(4.4)  one  need  only  evaluate 

F*  and  G*. 


covariance  (meters/sec) ~2 


*1 


Figure  4  -  Smoothing  error  for  velocity  as  a  function 


Figure  6  -  Pressure  at  1570  rads/sec 


Figure  7- Velocity  at  1570  reds/sec 
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